#!/usr/bin/env R

library(xts)
library(TSA)
library(vars)
library(sspir)
library(fGarch)
install.packages("sspir")

start_year <- function (data) {
  day <- index(data[1])
  year <- 1900 + as.POSIXlt(day)$year
  return(as.character(year))
}

end_year <- function (data) {
  day <- index(data[length(data)])
  year <- 1900 + as.POSIXlt(day)$year
  return(as.character(year))
}

both_years <- function (data) {
  return(paste(start_year(data), end_year(data), sep="::")) 
}

shift_dates <- function(dane, times, shift) {
  return(xts(dane[1:length(dane)], index(times[shift:length(times)])))
}

prepare_var <-function(dane, ...) {
  l <- list(...)
  l <- lapply(l, function(x) return(x[both_years(dane)]))
  return(l)
}

abs_min<-function(x, limit) {
  if(abs(x) < limit) {
    return(x)
  } else {
    return(limit * (x/abs(x)))
  }
}

limit_change <- function(data, limit) {
  if(is.na(limit)) {
    return(data)
  } else {
    return(as.xts(rollapply(na.omit(data), 1, function(x) return(abs_min(x, limit)))))
  }
}

do_diff <- function (dane, use_diff) {
  if(use_diff) {
    ret <- diff(dane)
    ret[1] <- 0
    return(ret)
  } else {
    return(dane)
  }
}

do_undiff <- function (result, dane, use_diff, max_change=NA) {
  if(use_diff) {
    return(na.omit(lag(dane + limit_change(result, max_change))))
  } else {
    return(na.omit(lag(result)))
  }
}

#calculate summ of absolute differences
calc_diff <- function(seq, data) {
  return(sum(abs(seq - data), na.rm=T))
}

calc_mean_diff <- function(seq, data) {
  return(mean(abs(seq - data), na.rm=TRUE))
}
#Simulates model y_t = y_{t-1} + e_t
random_walk <- function(dane) {
  return(lag(dane, 1))
}

#smooth data 
avg_data <- function(dane, ile) {
  ret <- as.vector(rollapply(dane, ile, mean))
  ret <- shift_dates(ret, dane, ile);
  return(as.xts(ret))
}

pred_ar_one <- function(dane, wstecz, typ="ols") {
  options(warn=-1)
  sim <- ar(dane, order.max=wstecz, method=typ, demean=FALSE)
  options(warn=1)
  gc()
  return((predict(sim, n.ahead=1)$pred)[1])
}

#predict ar model
pred_ar <- function(dane, wstecz, naraz, typ="ols", use_diff=TRUE, max_change=NA) {
  data <- do_diff(dane, use_diff)
  res <- as.vector(rollapply(data, naraz, function(x) return(pred_ar_one(x, wstecz))))
  res <-  shift_dates(res, dane, naraz)
  return(do_undiff(res, dane, use_diff, max_change))
}

pred_avg <- function(dane, ile, factor=1 ,use_diff=TRUE) {
  data <- do_diff(dane, use_diff)
  res <- avg_data(data, ile) * factor
  return(do_undiff(res, dane, use_diff)) 
}

#ar model on smoothed data
pred_ar_avg <- function(dane, wstecz, naraz, usrednij=5, typ="ols") {
  data <- avg_data(dane, usrednij)
  res <- pred_ar(data, wstecz, naraz, typ, use_diff=T)
  res <- ((res - random_walk(data)) / usrednij) + random_walk(dane)
  return(res)
}

pred_arma_one <- function(dane, wstecz_p, wstecz_q, t, use_diff) {
  options(warn=-1)
  sim <- tryCatch(arima(dane, order=c(wstecz_p, 1, wstecz_q), method=t), error = function(e) return(NaN), finally=1)
  if(is.na(sim)) {
    options(warn=1)
    print("Bad")
    #return random walk if there was an error
    if(use_diff) {
      return(0) 
    } else {
      return(dane[length(dane)]) 
    }
  } else {
    options(warn=1)
    return(predict(sim, n.ahead=1)$pred[1])
  }
}

pred_arma <- function(dane, wstecz_p, wstecz_q, naraz, typ="CSS-ML", use_diff=TRUE, max_change=NA) {
  data <- do_diff(dane, use_diff)
  res <- as.vector(rollapply(data, naraz, function(x) return(pred_arma_one(x, wstecz_p, wstecz_q, typ, use_diff))))
  res <- shift_dates(res, dane, naraz)
  return(do_undiff(res, dane, use_diff, max_change))
}

pred_var_one <-function(all, p, type, lag.max) {
  sim <- VAR(all, p=p, lag.max=lag.max, type=type);
  return((predict(sim, n.ahead=1)$fcst[[1]])[[1]])
}

pred_var <- function(dane, inne, naraz, p=1, type=c("both"), lag.max=NULL, use_diff=TRUE, max_change=NA) {
  data <- do_diff(dane, use_diff)
  other <- lapply(inne, function(x) return(do_diff(x, use_diff)))
  all <- data
  for (i in 1:length(other)) {
    all <- cbind(all, as.xts(other[[i]]))
  }
  j = 1;
  res <- vector(length=(length(data) - naraz + 1))
  for (i in 1:(length(data) - naraz + 1)) {
    res[j] <- pred_var_one(all[i:(i + naraz - 1)] , p, type, lag.max)
    j <- j + 1
  }
  res <- shift_dates(res, dane, naraz)
  return(do_undiff(res, naraz, use_diff, max_change))
}

pred_garch_one <- function (data, formula){
  sim <- garchFit(formula, data=data, trace=F)
  #gc()
  return(predict(sim)$meanForecast)
}

pred_garch <- function (dane, formula=~garch(1,1) , naraz, use_diff=T, max_change=NA) {
  data <- do_diff(dane, use_diff)
  res <- as.vector(rollapply(data, naraz, function(x) return(pred_garch_one(x, formula))))
  res <- shift_dates(res, dane, naraz)
  return(do_undiff(res, dane, use_diff, max_change))
}

pred_holt_one <- function (dane) {
  res <- HoltWinters(dane, gamma=F)
  return(predict(res))
}

pred_holt <- function (dane, naraz, use_diff=T) {
  data <- do_diff(dane, use_diff)
  res <- as.vector(rollapply(data, naraz, function(x) return(pred_holt_one(x))))
  res <- shift_dates(res, dane, naraz)
  return(do_undiff(res, dane, use_diff))
}

compare_rw <- function(prediction, data, year=2) {
  if(year==2) {
    ed <- end_year(data)
  } else if(year==1) {
    ed <- start_year(data)
  } else {
    ed <- year
  }
  ind <- index(prediction[ed])
  total_pred <- calc_diff(prediction[ind], data[ind])
  total_rw <- calc_diff(random_walk(data)[ind], data[ind])
  
  mean_pred <- calc_mean_diff(prediction[ind], data[ind])
  mean_rw <- calc_mean_diff(random_walk(data)[ind], data[ind])
  if (mean_rw >= mean_pred) {
    print("ZWYCIESTWO")
  }
  cat("    Całkowity błąd: \t predykcji=", total_pred, " \trandom walk=", total_rw, "\n")
  cat("    Błąd średni: \t predykcji=", mean_pred, " \trandom walk=", mean_rw, "\n")
  return(mean_pred)
}

check_cause <- function(dane, other, p=1) {
  razem <- cbind(dane, other[both_years(dane)])
  colnames(razem) <- c("dane", "other")
  v <- VAR(razem, p=p)
  print(causality(v, cause="other"))
  print(causality(v, cause="other", boot=T))
  print(causality(v, cause="other", vcov.=vcovHC(v)))
}

predict_zf <- function (data) {
  return(pred_ar(data, 1, 10))
}

predict_goog <- function(data) {
  return(exp(pred_ar(log(data), 1, 250)))
}

predict_cl <- function(data) {
  return(exp(pred_garch(log(data), ~garch(1, 1), 200)))
}

predict_eur <- function(data) {
  return(exp(pred_ar(log(data), 2, 500)))
}
